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IMPULSIVE ORBIT TRANSFER OPTIMIZATION 
BY AN ACCELERATED GRADIENT METHOD® 


By Ivan L. Johnson, Jr., Manned Spacecraft Center, NASA 


INTRODUCTION 


An analysis of impulsive orbit transfer optimization using an 
accelerated gradient method is presented. Several years ago in the 
Trajectory Optimization Section of the Mathematical Physics Branch, 
research began on parameter optimization methods especially for trajectory 
optimization. Since that time, several different methods have been tried 
and abandoned. Among those methods were the gradient projection method 
(ref. 1 and 2), the method of steepest descent using a penalty function 
(ref. l), various attempts at second order methods - all of which failed, 
and the conjugate gradient methods compared in reference 3. 

A quadratically convergent gradient method developed by W, C. Davidon 
(ref. U) was discovered. This method was refined slightly by Fletcher 
and Powell (ref. 5), and it was shown that it is probably the most 
powerful method known for the minimization of a general function of n 
variables . 

For the constrained minimum problem, it was found that if one minimized 
a quadratic type of penalty function (ref. l) for a given set of "large" 
penalty constants, an approximate solution can be obtained. If Davidon's 
method is used to minimize the penalty function, the inverse of the matrix 
of second partial derivatives of the penalty function at the minimum is 
obtained. If Newton's method is applied to the first-order necessary 
conditions for a constrained minimum, the exact solution for the constrained 
minimum problem can be obtained, provided the matrix of second partial 
derivatives of the augmented function can be calculated. Assuming that 
the approximate solution obtained by minimizing the penalty function using 
Davidon's method is "close" to the exact solution, the matrix of second 
partial derivatives needed for the solution by Newton's method can be 
approximated using the matrix from Davidon's method. This logic was used 
to develop a digital computer program known as the "Accelerated Gradient" 
program vhich is described in this paper. 

Presented at the Astrodynamics Conference held December 12, 13 and 
lL, 19^7, at the Manned Spacecraft Center, Houston, Texas. 
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Initially, the first type of orbit transfer models contained ellipses 
and circles only. Both the two-body and conic partials equations were 
good for ellipses and circles only. In order to solve problems of the 
more complex form involving hyperbolas as well as ellipses and circles, the 
universal variable formulation by S. Pines (ref. 6) and A. K. Nakashima 
(ref. 7) is now used. This formulation is also described briefly in 
this paper. 

The formulation, programming, and accumulation of data for a three- 
impulse orbit transfer problem was accomplished mainly by W. C. Bean of 
the Mathematical Physics Branch, NASA, Houston, Texas, using the 
Accelerated Gradient program with the universal variable, two-body formu- 
lation. The optimal solution to thi* problem is presented in this paper 
as an example. In the problem, both an intermediate trajectory constraint 
and an inequality constraint are featured. In a report to be published, a 
further description and more data on this example problem will be presented 
by Mt". Bean. 


i 

STATEMENT OF IMPULSIVE ORBIT 
TRANSFER OPTIMIZATION PROBLEMS 


With impulses approximating the rocket burns and two-body motion 
(conics) approximating the unpowered flight, orbit transfer models of 
various missions can be constructed for the purpose of parameter 
optimization. 

The performance index for orbit transfer optimization problems is 
often the "characteristic velocity," the sum of the magnitudes of the 
impulsive-velocity vectors. The constraints, inequality and equality 
types, are usually specified on the terminal orbit, but constraints may 
be placed on intermediate trajectories as well. 
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The generalized eccentric anomalies, which are the independent 
variables governing the lengths of the coasting arcs, and the components 
of the impulsive-velocity vectors are the parameters of the optimization 
problem. 

The orbit transfer optimization problem then is to find the param- 
eters minimizing the performance index subject to the system of trajectory 
constraints . 


A DIGITAL COMPUTER PROGRAM 
FOR PARAMETER OPTIMIZATION 

The Accelerated Gradient program, a digital computer program, using 
the method of reference 8, is used for the impulsive orbit transfer 
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optimization. The program numerically solves the following general 
optimization problem. 

• Find that set of parameters , ®2***** a n w ^ich locally minimizes 
the general function 




f = f(a^» oigi ...» a^) 

subject to the system of m < n nonlinear constraints 


( 1 ) 




a 2* •••» * 0 
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S 2^ a l* a 2* **** a n^ = 0 
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Included in this system are inequality constraints. 

if 




h,(a , a 0 , ...» a h <_ 0 
J 1 2 n ~ 


(3) 

• • 


is needed, then with the help of the function S(hj), 
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S(h ) =0 if h < 0 
J J 

9 


(U) 
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S(h ) =1 if h >0 
J o 





inequality constraints can be redefined as equality constraints 



gj = hjS(hj) = 0 where J <_ 

: m * 

(5) 



The Accelerated Gradient program consists of two phases. 

The first 



phase treats the optimization problem in an approximate form, 
unconstrained function (penalty function) 

The 



m 
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f = f + iLv j 2 


(6) 
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1 8 formed and minimized for given "large" penalty constants k > 0. 
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Should solutions of both the approximation problem and the original 
problem exist, the former approaches the latter as each • (ref. 1). 


The method used for minimizing f was developed by William C. Davidon 
(ref. U, 5, and 9). Davidon* s method is a quadratically convergent 
gradient method for the minimization of a general function of n variables. 
Beginning with a starting point (a^, a 0 , ...» a^), the gradient vector 


and the function 


are calculated. Using the formula 


Aa = -yHf 


a change in the vector a of parameters a^, a^, . . . , is made. The 

matrix H is symmetric and positive definite (it must be chosen as such 
initially), and the scalar y is a positive step-size parameter. The 
one-dimensional minimum of f(o + Ao) versus y is obtained (the method 
of reference 10 is used in the computer program), the gradient vector 
? and the function f is recalculated at the new a, and H is updated ac- 


cording to the formula 


H + AH = H + 


AaAa' 


HA? A? T « 
at a 


T - 

Aa Af 


Af T HAf 
a a 


The procedure is then repeated using the new values of s f a » and H until 

the set of parameters (a^, a 0 , ...» a^) is found which locally minimizes f 

for a given set of penalty constants k > 0. As a bonus, upon covergence 

J 


to the minimum of f, the H matrix approaches the inverse of the matrix 
of second partial derivatives of f with respect to the parameters 


a at the minimum. That is, 
J 


\ u 'a*a 
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In order for the solution to the approximation problem to approach 
that of the original problem, the penalty constants k. may be enlarged 
according to the relationship ^ 




( 10 ) 


and f is minimized again. 




fefore repeating the minimization of f with larger values of the 
the initial H matrix will be calculated using the converged one 


from the previous minimization and updating it according to the changes 
in the k^ . The initial H matrix is found by iterating m times using 

the formula 


H «- H - 



Ak. 


1 + A* l6ici Hg l0 



( 11 ) 


V . 


Mb 


SL- 

Ml 




where the symbol «- denotes "is replaced by" and g^ is the gradient vector 

of the i-th constraint at the previous minimum. The function f is then 
minimized again. 

The second phase of the Accelerated Gradient program consists of taking 
the solution found by the first phase and refining it using an algorithm 
obtained by applying Newton’s method. 

The algorithm is obtained as follows. The first order necessary 
conditions for the constrained minimum problem is given by the system 
of equations 


f + g X 
a a 


g = 0 


I 

i ‘ 


( 12 ) 


■» 






where g^ is an n * m matrix whose columns are the gradient vectors of 

the constraints, X is an m-dimensional vector of multipliers and g is 
an m-dimensional vector of the constraints. Application of Newton's 
method to the system of equations (12) yields 
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(f ♦ x T g) g 


da 


f ♦ g x 

aa a 




a a 




■ « 


g„ T 0 


dX 


g 

a 
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(13) 


Using the converged H matrix* the matrix of second partial derivatives 
T 

(f + X g) Qa can be approximated by 


H-> - g 0 Kg‘. (1») 

where K is an m * m diagonal matrix, the diagonal elements being the 
penalty constants k,. In terms of known quantities, the expression 

v 

for da can be obtained from equation (13) as 

m .1 

(15) 


do - -[h - Hg 0 ( 6 X>'V H ] f 0 - H 8 0 (g o TK8 « ) ' lg - 

Re, ..ated application of equation ( 15 ) starting with the "near” solution 
obtained by minimizing f using Davidon's method yields a solution 
to the original problem. Convergence is obtained when the da^'s are con- 
sidered small enough. » 


TWO-BODY MOTION USING A UNIVERSAL VARIABLE FORMULATION 

For the coast phases of the orbit transfer models, the two-body 
equations are presented in a form using a universal variable formulation. 
The advantage of this type of formulation is that only one set of two-body 
equations is needed to describe all of the different conics. 


The generalized eccentric anomaly 6 is defined by 

yB 2 = e 2 


where 


. 1 = _2 _ _o_ 

Y " a ~ r ~ u 
o 


(16) 


(IT) 


and r Q is the magnitude of the initial position vector. R q , v q is the 
magnitude of the initial velocity vector, V q , and \i is the gravitational 




constant for the central body. For elliptical orbits 

e • E - e . 


UO) 


where E is the eccentric anomaly. For hyperbolic trajectories 0 is 
imaginary. 

The equations of motion in a central force field are given by 


-px 

x t * , i * 1, 2, j 

r 3 


(19) 


where r is the magnitude of the position vector R, x.^ (i * 1, 2, 3) 

are the components of R, and (i * 1, 2, 3/ are the components of V. 

With the initial position vector, R , the initial velocity vector, V , 

o o 

and the initial time, t Q given, integrals of the equations of motion 
(19) with 6 as the independent variable are given by 


X 1 * fx oi + gx ol 'i 

a . a a ) » 

X 1 * fx ol + gx oi J 


1 - 1, 2, 3 


The coefficients f, g, f, and g are given by 


(20) 


f = 1 - — , 
r 

o 


— r 0 + -Sg J] , 

G (_ ° 1 G 2 J 


• 47 

f = - G, , 
rr 1 
o 


( 21 ) 


(22) 


(23) 




(24) 




I 


where 


d - R • V 
o o o 


r ■ r G < ' 7 "G ♦ • 

o o /y 1 2 


(26) 


The universal variables G p (p * 0, 1» 2, ...) are defined by 


v 1 e 4 

° P * 6 FT ' ( P VsTT + ( P VTJT 


• • 


(2T) 


Calculating the G p with the highest valued subscript, say G p and 

G . , by the series in equation (27), the remaining G ’ s 
P-1 P 

(G G G ) can be calculated using the recursive equation 

p-2 p-3 o 


J p-2 * (p - 2)! Y °p' 


(28) 


The generalized Kepler’s equation is given by 


t = t + — G 0 + r G. + -^ G-. 
0 /y 3 0 1 ft 2 


(29) 


With the initial state, R , V q , and t Q , given, the state of an orbit, 

P, V, and t, for any specified 6 is given by equations (20) and (29). 
Determining the state in this manner gets rid of iterations involving 
Kepler's equation. 

A detailed account of the equations presented in this section is 
presented in references 6, 7, 11, 12, and 13. 


TRAJECTORY DERIVATIVES 

In using a gradient method for constrained minimization, the gradient 
vectors of both the function and the constraints are calculated. For 
orbit transfer optimization problems, the partial derivatives of both 


' ,v 
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the fYnction and constraints with respect to the components of the impul- 
cive velocity vectors and the generalized eccentric anomalies must be 
calculated. For the gradient vectors of the constraints, the partial 
derivatives of the state at one time with respect to the state at an 
earlier time and the derivatives of the state with respect to the 
generalized eccentric anomaly are needed. 


From equations (l6) through (29), the state transition matrix Q 
c f partial derivatives of the state at one time with respect to the 
£tu.te at an earlier time can be derived in closed form. 


Q5 


m 

(•Sr) 


Jiyj(%) 
&)!(*) 


(t)K^) 


(30) 


The elements of the 3*3 submatrices (3 x^/3x q j), (3x^/3x q j), (3x^/3x^j), 


oj 


oj 


and (3x./3x ), are given by 

1 O J 






V 


L-Jki V< 

ij r r • 


X . X 4 1 

-2J.0 +-2J- G, 

r o /—■ 1 

i o * w i 


Vi I /, . °3 - BG 2\ 

1 ^ / 


+ ^ ° 


♦ (2G|. - 6G,) x 


[ ol * '^U * Bu 3 ' x oi * 


(33) 


g6 ij + 


k - *oi> 


• • 
^2-21-0,* , ♦ -2i 
/T 1 °J ur 

mJ V* 


(G 3 ’ aG 2 } ^ X oi 
o 


+ (aoj. - » 3 > Vi 


(3*0 


The elements of the 3 * 1 submatrices (3x. /3t ) and (3x/3: ) are 

1 o o 

given by 


3x^ 3x^ 


3t 3t 
o o 


= 0, 


(35) 


The elements of the 1 * 3 submatrices (3t/3x .) and (3t/3x .) are given by 

oj oj 


3t .. *£. (i • ) , 

3 x oJ ^3 P 1 « J J 


(36) 


= - — (p x - G x ), 
W 1 0.1 2 oj 


(37) 


= 1. 


(38) 


In the preceding equations, 6^ is the Kronecker delta. 


1, when i = J \ 

0, when i # J j 


( 39 ' 


\ 
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and 


r 

p 2 " " [ r o 6G l " ^ 


(3G 5 - 6G U ) + r Q (G 3 - 6G 2 ) + (20^ - 6G 3 ) 




(G 3 - 6G 2 ) - (2G U - 6G 3 ) 


• (U0) 


(4l) 


The derivatives of the state with respect to the generalized eccentric 
anomaly are given by 


• < 



(42) 


and 



dt _ r 

dB ~ /u ‘ 


(43) 


(44) 


AN OPTIMAL THREE-IMPULSE TRANSFER 

The orbit transfer optimization method presented was applied to a 
three-impulse orbit transfer featuring an intermediate constraint. 

The orbit transfer that was optimized may be described as follows. 
Starting with the initial state, R q , V q , and t Q , with respect 

to the earth as the central body, an initial impulsive velocity vector 
AV & is applied, followed by a coast, whose length is governed by 6 Q , 

on an ellipse whose period P is governed by the constraint 

g 1 = P - P = 0. (45) 

Another impulse, AV 1 » is applied followed by a coast whose length is 
governed by 6 . A third impulse, AV 2 » is applied to transfer onto a 
hyperbola, partially specified by the constraints 
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g 2 = D - D = 0 


g_ = v - v =0 
°3 00 00 


g, = r - r <0 

*1* p p - 


) ’ 


(46) 


where D and v^ are the angle of declination and magnitude of the 

velocity vector at infinity, respectively. The quantity r^ denotes 

the magnitude of the radius vector at perigee, and (~) denotes given or 
specified quantities. 

The parameters selected for optimization are 

^2 
Az„ 


a = 


Ax 


1 


Ay-, 


Az. 


Ax 


Ay, 


Az 


(47) 


The performance is measured by 

f(«) “ |av q | + |av x | + |av 2 |. 


(48) 


The function f is minimized subject to the constraining equations (1+5) 
and (46) where 


( 49 ) 


v = 37 002.647 fps, 
00 * 


r = 35U3.93U n. mi. (100-n. mi. radial altitude). (50) 

P 

The initial state conditions are those of a circular orbit with a radial 

altitude of 262.0 n, mi. The vectors R and V are given by 

o o 


3705. 93U n. mi. 


(51) 


V q = 25 002. 647 fps , 


(52) 


and t =0.0 seconds, 
o 

Figure 1 is a pictoral display of the case where D 
optimum transfer is a double Hohmann transfer. 


= 0. The 
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Figure 1.- Optimized three impulse orbit transfer. 
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Figure 2 is a graph of the performance f versus the angle of declina- 
tion D for different values of the period P of the first transfer 
ellipse. Observing the graph, the performance gets better as the period 
of the second ellipse increases and gets worse as the angle of declination 
of the velocity vector at infinity for the hyperbola increases. 


NASA 






IMPULSE 

VELOCITY 



Figure 2.- Impulse velocity versus angle of declination 

for given periods . 


One interesting result worth mentioning is that for the cases where 
D ^ 0, the impulses did not occur at the apsidal points of the ellipses 
or perigee of the hyperbola. Also, the impulses were not along the same 
line as the velocity vectors. For D = 0 the impulses did occur at the 
apsidal points and the impulses were in the same direction as the velocity 
vectors . 


CONCLUDING REMARKS 


For the orbit transfer models considered in this paper, the perfor- 
mance indices were defined as the sums of the magnitudes of the impulse 
velocity vectors. This seems to be the more common type of performance 
index; however, it could be defined as some other quantity. Also, 
besides the parameters mentioned, other variables could be chosen as 
parameters as well. 


The reason for using a universal variable formulation for the solution 
to the two-body equations is that the single set of equations defining 
the solution is good for all the different types of trajectories (conics). 


L. 
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Specifying the state of a trajectory with seven dependent variables (the 
three components of both the position vector and velocity vector, and time) 
rather than the usual six gets rid of the usual iterations involving 
Kepler's equation. A good account of the numerical evaluation of the 

infinite series for G is given in reference 7 • 

P 

The Accelerated Gradient program is not only used for orbit transfer 
optimization, but because it is written in a general form, it can be 
applied to optimization problems of a different nature. For example, the 
program is currently being used for the optimization of an electronics 
system and a propulsion system. 


* 
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